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The presence of gaseous giant planets whose 
orbits he in extreme proximity to their 
host stars ("hot Jupiters"), can largely be 
accounted for by planetary migration, as- 
sociated with viscous evolution of proto- 
planetary nebulae^. Recently, observations 
of the Rossiter-McLaughlin effect^ during 
planetary transits have revealed that a con- 
siderable fraction of detected hot Jupiters 
reside on orbits that are misaligned with 
respect to the spin-axes of their host stars^. 
This observational fact has cast significant 
doubts on the importance of disk-driven mi- 
gration as a mechanism for production of 
hot Jupiters, thereby reestablishing the ori- 
gins of close-in planetary orbits as an open 
question. Here we show that misaligned 
orbits can be a natural consequence of disk 
migration. Our argument rests on an en- 
hanced abundance of binary stellar compan- 
ions in star formation environments, whose 
orbital plane is uncorrelated with the spin 
axes of the individual stars^'^'^. We analyze 
the dynamical evolution of idealized proto- 
planetary disks under perturbations from 
massive distant bodies and demonstrate 
that the resulting gravitational torques act 
to misalign the orbital planes of the disks 
relative to the spin poles of their host stars. 
As a result, we predict that in the absence 
of strong disk-host star angular momen- 
tum coupling or sufficient dissipation that 
acts to realign the stellar spin axis and the 
planetary orbits, the fraction of planetary 
systems (including systems of hot Neptunes 
and Super-Earths), whose angular momen- 
tum vectors are misaligned with respect to 
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Fig. 1. — Geometrical setup of the problem. This fig- 
ure depicts a schematic representation of the produc- 
tion of misaligned close-in planets through disk-driven 
migration in binary systems. The adiabatic response 
of a self-gravitating disk to long-term perturbations 
by a stellar companion is the recession of its ascend- 
ing node, as defined by the orbital plane of the stellar 
companion. The recession of the disk's angular mo- 
mentum (AM) vector about the stellar binary orbital 
AM vector appears to be an excitation of misalign- 
ment between the stellar spin- axis and the disk, -0, in 
the star's reference frame. 



their host-stars should be commensurate 
with the rate of primordial stellar multi- 
plicity. 

The obliquities (angles between the planetary 
orbits and the stellar spins) of detected plane- 
tary orbits range from almost perfectly aligned 
prograde to almost perfectly aligned retrograde 
systems^. Previously, the misalignment between 
planetary orbits and stellar spin axes had been at- 
tributed to post-nebular multi-body interactions. 
Most notably, Kozai cycles with tidal friction^'^'^^, 
planet-planet scattering^ ^'^^, and chaotic secular 



excursions'^ have been invoked as a means of pro- 
ducing misaligned planets. These mechanisms are 
likely responsible for a few specific examples (e.g. 
the extreme eccentricity of HD80606b is almost 
certainly due to Kozai resonance with the stellar 
companion HD80607^), however it is unlikely that 
they can explain misaligned hot Jupiters as a pop- 
ulation. For example, the Kozai mechanism can be 
stified by forced apsidal precession in multi-planet 
systems'^. Likewise, within the context of planet- 
planet scattering and secular chaos, the allowed 
parameter range is limited, since the production of 
close-in orbits requires the timescale for tidal cap- 
ture to be considerably shorter than that for eccen- 
tricity growth'^, while demanding the associated 
tidal heating to be sufficiently small to not over- 
inffate the planet beyond its Roche- lobe'^. Ad- 
ditionally, the observed presence of mean-motion 
resonances among giant planets on wide orbits 
(which rely on smooth, convergent migration to 
congregate'^), provides further motivation for the 
development of a unified model for disk migration 
that is capable of producing misaligned orbits. 

The dynamics of self-gravitating proto-planetary 
disks under external perturbations can be ex- 
tremely complex, making precise quantitative 
modeling computationally unfeasible. Conse- 
quently, here we concentrate on characterization 
of the qualitative physical behavior of the system 
and utilize classical perturbation methods to ob- 
tain a solution. In the spirit of secular theory'^, 
we model the proto-planetary disk as a series of 
initially planar, circular, concentric massive wires 
that interact gravitationally. Our model is based 
on the Gaussian averaging method'^' '^ and the 
gravitational potential is softened in order to par- 
tially account for the discrete representation of 
the disk. The effects of dissipative ffuid forces 
within the disk are neglected. The perturbing 
body is also modeled as a massive ring, but is ec- 
centric {e' = 0.5) and inclined with respect to the 
disk by an inclination i' . A detailed description 
of the model and its inherent assumptions is pre- 
sented in ref.[19] while the particularities of our 
implementation are stated in the Supplementary 
Information (SI). 

A self-gravitating disk will preserve an un- 
twisted structure and act as a rigid body, pro- 
vided that the characteristic timescale of the ex- 
ternal perturbation greatly exceeds that of the 
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Fig. 2. — Excitation of disk-star misalignment. The 
time evolution of the misalignment angle, ip is shown. 
The considered rridisk — 10~^Msun — lOMjup nebula 
has a surface density profile of the form S oc r~', 
extends between ain — 1 AU and aout — 50 AU, 
and orbits a Mstar = IMsun host star. The plot- 
ted curves represent the dynamical states of the disk 
annuli, where every other ring is plotted. The gray 
curves depict the inner-most annuli. The close prox- 
imity of the rings to to each-other demonstrates that 
the disk remains locally un- warped and acts as a rigid 
body to a good approximation. The extent of rigid- 
ity is largely controlled by the disk mass, with heavier 
disks remaining closer to their mid-plane. The host 
star is assumed to be a slow rotator and its spin pole is 
held fixed for simplicity. The perturbing M' = IMsun 
binary companion lies a' = 500AU away, has an ec- 
centricity of e = 0.5, and is inclined with respect to 
the disk by z = 45° (red) and by by z = 70° (blue). 
Throughout the duration of the integration, the an- 
nuli of the disk never attain significant eccentricity 
(e < 0.1 - see SI for an in-depth discussion). This 
calculation was performed using a conservative soft- 
ened Gaussian averaging model (with Airings = 31 - 
see SI), and thus contains no restrictions on the sec- 
ular dynamics of the system, but ignores dissipative 
forces of the gas. The nodal recession periods charac- 
teristic of this setup are Tdisk — 0.9Myr {i = 45°) and 
Tdisk ^ l.SMyr (i = 70°). The results in this figure 
can be translated to other system parameters by not- 
ing that the the maximum misalignment attainable by 
the disk is roughly twice the disk/binary orbit incli- 
nation and that the recession period scales approxi- 
mately as Tdisk oc a^^/(M^ cos i^(l + 3e^^/2)) (see also 
Figure 3). 



disk's self-interaction^^. Mathematically, this 
amounts to a statement of adiabatic invariance 
of the phase-space area occupied by a single sec- 
ular cycle within the disk^^. If this condition is 
satisfied, the external perturber's sole effect is to 
induce a recession (i.e. drift) of the ascending 
node of the disk, as defined by the plane of the 
stellar orbit. The embedded planetary orbit will 
also adiabatically follow the disk. 

In the reference frame of the host star, the nodal 
recession of the disk will appear as a cyclic excita- 
tion of inclination between the disk and the stel- 
lar spin axis (see Figure 1), provided that the host 
star's angular momentum vector does not adiabat- 
ically trail the disk. For this to hold true, the char- 
acteristic interaction timescale between the disk 
and the stellar spin-axis, Tstar must exceed the 
disk's nodal recession timescale, Tdisk^ by a con- 
siderable amount (i.e. angular momentum cou- 
pling between the disk and the host star must be 
non- adiabatic)^. The former can be estimated by 
modeling the stellar rotational bulge as an iner- 
tially equivalent orbiting ring, effectively reduc- 
ing the characteristic interaction timescale to the 
forced nodal recession period of the ring . 

Observations suggest that rotational periods of 
T Tauri stars, whose masses exceed M > O.25M0, 
form a bimodal distribution where fast and slow 
rotators are centered around ~ 2 days and ~ 8 
days respectively, with a preference for slow ro- 
tation at higher masses^^. Thus, for typical pre- 
main-sequence stars, we obtain Tstar ^ 10 Myr 
and Tstar ^0.3 Myr for slow and fast rotators re- 
spectively (see SI for details). As will be shown 
below, this suggests that the adiabatic trailing of 
the stellar spin axes will only prevent excitation 
of mutual misalignment for fast rotators. Further- 
more, the disk-star angular momentum transfer 
will likely be unimportant in mature disks because 
of low accretion rates^^. 

In addition to avoiding the adiabatic trailing 
of the host stars, the prominence of the mecha- 
nism described here is determined by the abun- 
dance and longevity of wide stellar binary sys- 
tems in star formation environments, since the 
misalignment angle, ?/^, becomes frozen in when 
the binary companion is stripped away or when 
the proto-planetary disk dissipates (nevertheless, 
a long-lived binary companion can act to misalign 
a mature planetary system with its host star^^). 



While it is tempting to estimate the frequency 
with which significant misalignments are attained 
by this mechanism via population-synthesis, the 
enormous range and vast observational uncertain- 
ties in the input parameters would render such 
a calculation of little practical use. In particu- 
lar, although observations of the Taurus- Auriga 
star-forming region^^ suggest that the orbital dis- 
tribution of young solar-type binaries is roughly 
log- flat with an overall binary fraction of ~ 40%, 
it is noteworthy that the process of wide binary 
formation also appears to exhibit environmental 
dependence^. Simultaneously, the rate at which 
wide binaries get disrupted in birth clusters de- 
pends sensitively on the local densities within the 
clusters^, which remain observationally elusive, as 
the majority of stars are born in aggregates that 
dissolve quickly (on timescales of a ~ few x 10 
Myr or less)^^. Still, the above conditions likely 
imply that the timescale for excitation of signif- 
icant misalignment should be considerably less 
than -- lOMyr. 

In systems where self-gravity is strong enough 
to maintain the effective rigidity of the disk, 
fast circulation of the disk's argument of peri- 
helion also ensures adiabatic eccentricity dynam- 
ics. Specifically, this means that the disk will 
not develop significant eccentricity, as the Kozai 
resonance within the individual annuli will be 
suppressed^^ (see SI). Recall also that here, we 
are ignoring dissipative effects that would gener- 
ally act to circularize the disk and maintain its 
rigidity. Assuming that the angular momentum of 
the stellar binary greatly exceeds that of the proto- 
planetary disk, the maximal inclination that can 
be excited between the stellar rotation axis and 
the disk's orbital plane is approximately equal to 
twice the inclination of the proto-planetary disk 
with respect to the binary orbit: tl^max — 2^^ As 
a result, retrograde planetary orbits can be nat- 
urally achieved, provided that the inclination of 
the stellar companion exceeds 45 degrees. Figure 
2 shows the time-evolution of two such examples. 

While the Gaussian averaging model employed 
above yields a rigorous representation of the secu- 
lar evolution of the system, it is also computation- 
ally expensive. Fortunately, similar results can be 
obtained with a modified (arbitrary i^ e') Laplace- 
Lagrange analytical theory ^^ (see SI), providing an 
avenue for efficient mapping of parameter space. 
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Fig. 3. — Timescales for excitation of spin-orbit mis- 
alignment. The characteristic nodal recession period 
of a rridisk — 1^~^ Mstar disk, as a function of binary 
mass and orbital properties is shown. The considered 
disk has an outer edge at aout — 50 AU. Increasing 
dout will result in linear increase of the precession fre- 
quency. Note that the period is expressed as a scaling- 
law in the mass of the host star. In the region of pa- 
rameter space where the precession period greatly ex- 
ceeds the disk lifetime, only small misalignment angles 
between the disk and the host star can be excited. In 
principle, however, if the stellar companion does not 
get stripped away, the ascending node of the invari- 
able plane of the formed planetary system can also 
recess. However, the degree to which this can affect a 
planet on a close-in orbit is sensitive to the particular 
architecture of the system. In the region of the pa- 
rameter space where dynamics seizes to be adiabatic, 
misalignment is certainly attainable, but more quan- 
titatively precise (magnetohydrodynamical) modeling 
is required for its characterization. Finally, on the 
extreme high-mass/small orbital separation end of pa- 
rameter space, one could envision a scenario where a 
newly formed disk becomes severely twisted and even- 
tually gets disrupted as a result of strong external per- 
turbations. Driven by viscous dissipation, however, 
such a structure would likely re-collapse into a new 
protoplanetary disk, whose orbital plane will be close 
to the Laplace plane of the stellar binary orbit. 



We have quantified the precession timescale for 
a range of stellar companion masses as well as 
binary separations. The results are summarized 
in Figure 3. This calculation highlights an ef- 
fective equivalence between distant massive per- 
turbers and lower mass perturbers with smaller 
semi- major axes, since the recession period of 
the disk scales approximately as Tdisk ^ a'^ /M' . 
Thus, the precessional effect of a M' = IMq star 
orbiting at a = 10^ AU is equivalent to the pre- 
cessional effect that arises from the protoplane- 
tary disk and its host star orbiting a ~ 10^ M© 
star cluster at a = 0.25 pc. It is further note- 
worthy that bound companions are not necessarily 
required for production of oblate disks, as impul- 
sive perturbations from passing stars in the birth 
cluster will cause the inclination of the disk to ex- 
ecute a random-walk and in some cases can ex- 
cite significant misalignment^^. Collectively, this 
explanation places the ~ 7 degree misalignment 
between the Sun's spin-axis and solar system's in- 
variable plane into a more general, extrasolar con- 
text. Although, the process of planet formation in 
perturbed, warped disks certainly deserves further 
study. 

Given the diverse nature of the environments 
in which planetary systems may form, one would 
expect a wide range of characteristic precession 
timescales for proto-planetary disks. Although 
this does not necessarily imply an isotropic distri- 
bution of spin-orbit angles, it is quite possible that 
hot Jupiters which emerged from protoplanetary 
disks in multi- stellar systems, already resided on 
misaligned orbits at the time of nebular dispersion. 
If this is true, small spin-orbit angles must be in 
large part either a result of adiabatic trailing of the 
host stars or dissipative re-alignment of the sys- 
tem. Such a scenario appears to be supported by 
observations: hot Jupiters orbiting hot, massive 
stars tend to be misaligned while hot Jupiters or- 
biting less massive cooler stars tend to have small 
spin-orbit angles^. This transition has been at- 
tributed to tidal re-alignment of cooler stars due 
to the increased size of their convective zones^^ 
and thus appears to be in good agreement with 
our model. 

The consistency arguments presented above 
indicate that disk-driven migration in binary 
systems is a favorable origin of misaligned hot 
Jupiters. We can test this model as follows. Al- 



though multiple explanations exist for the origins 
of hot Jupiters, short-period multi-planet systems, 
which are typically less massive, have almost cer- 
tainly undergone disk-driven migration. Theo- 
retically, disk-driven migration tends to maintain 
coplanarity and near-resonances in the planetary 
systems^^, both of which have been spectacu- 
larly confirmed by the Kepler mission^^. Within 
the context of the model proposed here, proto- 
planetary disks become misaligned with the spin 
axes of their host stars irrespectively of the masses 
and orbital radii of the newly formed planets. 
We predict that future observations of Rossiter- 
McLaughlin effect will reveal that systems of close- 
in coplanar sub-giant planets can be misaligned 
with respect to the spin axes of their host stars 
in the absence of significant dissipative processes. 
Similarly, the presence of distant, nearly circu- 
lar resonant planet pairs in systems that host 
hot Jupiters on oblique orbits, would also point 
to disk-torquing as the likely mechanism for the 
origin of spin-orbit misalignment. High-precision 
radial velocity monitoring of transiting systems 
should directly test this prediction in the near fu- 
ture. 
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Supplementary Information 

In this work, our primary focus lies in understanding the 
long-term gravitational interactions between protoplane- 
tary disks and distant inclined external perturbers that are 
gravitationally bound to the host-stars. We model the disk 
as a series of prograde-orbiting concentric, massive wires 
that are initially coplanar and nearly circular (e ~ 10~^). 
We average over the Keplerian mean motion of the exter- 
nal perturber and also model it as a massive wire. In the 
example shown in the main text, the perturber is taken to 
have an eccentricity of e^ = 0.5, and is inclined with respect 
to the disk by an inclination i' = 45 deg and i' = 70 deg. 

We work at two levels of accuracy: first we utilize an un- 
restricted (arbitrary e, i) numerical model based on Gaus- 
sian orbital averaging method. Second, we obtain quantita- 
tively similar results with an analytical small-angle approx- 
imation to the Gaussian model, i.e. the Laplace-Lagrange 
secular theory. 

1. Gaussian N-Ring Model 

Qualitatively, Gauss's averaging method is rather in- 
tuitive: calculation of the phase-averaged interactions be- 
tween two non-resonant orbits is equivalent to treating the 
two orbits as massive wires, where the line-density is in- 
versely proportional to orbital velocity, and computing the 
forces they exert on each-other. Within the context of this 
approximation, the semi-major axes of the rings become 
constants of motion (see Ch.7 of ref. [18]). This method 
allows for a rigorous treatment of a conservative represen- 
tation of the protoplanetary disk as well as the distant per- 
turber provided that the number of such wires is large. 
In practice, the full Gaussian averaging procedure can be 
computationally expensive, severely limiting the number of 
rings in the model. Furthermore, in a conventional calcula- 
tion, ring-ring crossing produces a discontinuity in the force 
calculation that requires high resolution in the integration. 

Fortunately, a recent extension of Gauss's method has 
allowed for significant advances in the implementation of 
the algorithm described above-*^^. In particular, ref. [19] 
extended the averaging method to softened gravitational 
interactions. By extension of gravitational softening of a 
point-mass, which is equivalent to dispersing the mass over 
a Plummer sphere*^-*^, softening the gravitational potential 
in an N-ring code allows for an approximate representa- 
tion of a continuous disk with a limited number of dis- 
crete wires. Furthermore, the analytic averaging approach 
demonstrated in ref. [19] results in a tremendous speed-up 
of the integration. 

The model employed in this work borrows directly from 
ref. [19]. In interest of conciseness, we shall not restate 
the details of the the softened Gaussian averaging method 
here and restrict ourselves to the particularities of our im- 
plementation. We model the disk as a system of 30 equal- 
mass rings, distributed in equal semi-major axis intervals 
(i.e. E oc r"-*^ surfrace density profile) between ain = lAU 
and Gout = 50 AU, comprising rudisk = 10~^ Mstar- We 
choose a softening length of e = 0.79AU (approximately 
half the separation between the rings). In the numerical 
averaging procedure, the rings are broken up into at least 
128 sectors each. The equations of motion are integrated 
using the fourth-order Runge-Kutta method"^^ with a con- 
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Fig. 4. — Eccentricities of the disk annuli as functions of 
time corresponding to the Gaussian secular solutions pre- 
sented in Figure 2. Although the disk is perturbed by a 
massive (M' = Mgtar) binary companion with e' = 0.5, 
the predominantly adiabatic nature of the perturbations 
ensures that significant eccentricity is never excited and the 
Kozai resonance is erased. As in Figure 2, the red curves 
corresponds to the solution with i' = 45 deg where as the 
blue curves correspond to the solution with i' = 70 deg. 



stant timestep of r = 100 years. Increasing the number 
of rings by a factor of two while decreasing the softening 
length or decreasing ain did not modify the results signifi- 
cantly. 

The primary effect of disk self-gravity is to offset the ra- 
dial frequency from the orbital mean motion, and thereby 
give rise to comparatively fast apsidal precession^^. This 
effect is of great importance to planetary formation as 
it acts to effectively decouple the eccentricity dynamics 
of the disk from that of the perturbing body. Specifi- 
cally, this is accomplished as follows. To leading order, 
secular excitation of orbital eccentricity is governed by a 
term in the Hamiltonian that contains a harmonic of the 
form cos('ci7 — -uu'), where vu is the longitude of perihelion 
and as before, the primed quantity refers to the external 
perturber-*^^. An application of Lagrange's planetary equa- 
tions yields de/dt oc sin(ti7 — m'). It is easy to show that 
significant excitation of the eccentricity can only be ac- 
complished if (m — m') is a slowly varying quantity (since 
the integrated effect of this harmonic scales inversely with 
d{vu — vu')/dt). Otherwise, the above mentioned harmonic 
becomes a quickly varying term, and has a small effect on 
the averaged evolution of the orbits (recall that the sec- 
ular approximation itself is motivated by filtering out all 
quickly- varying effects from the Hamiltonian). In direct 
analogy, other higher order terms the Hamiltonian con- 
taining longitudes of perihelion will also only act to intro- 
duce low-amplitude, high-frequency noise in the eccentric- 
ity solution. In other words, because the self-interaction 
timescale within the disk is much shorter than the exter- 
nal perturbation timescale, the eccentricity dynamics of the 
disk predominantly exhibit adiabatic effects^^. 

Aside from direct secular excitation of eccentricities, an- 
other seemingly detrimental effect which is erased by disk 
self-gravity is the Kozai resonance"^^ . An important prop- 
erty of the Kozai resonance is the libration of the argu- 
ment of perihelion. If the argument of perihelion of the 
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Fig. 5. — Time evolution of the stellar-spin/disk orbit 
normal misalignment angle. The setup is identical to that 
considered in the main text, with the exception of the bi- 
nary companion's eccentricity (here taken to be e' = 0.1) 
and inclination (here taken to be i' = 45° (red) and by 
i' = 75° (blue)). This calculation was performed using 
a conservative softened Gaussian averaging model (with 
brings = 31), and thus contains no restrictions on the secu- 
lar dynamics of the system, but ignores dissipative forces of 
the gas. The nodal recession periods characteristic of this 
setup are Tdisk - 1.4Myr (i' = 45°) and Tdisk - 4.6Myr 
(i' = 70°). Throughout the duration of the integration, 
the annuli of the disk never attain significant eccentricity 
(e < 0.01). 



disk circulates quickly (as it does in the calculations pre- 
sented here), the Kozai effect seizes to be a resonance and 
instead simply becomes another quickly-varying term in 
the Hamiltonian^'-'. These arguments are directly tested 
in the Gaussian N-ring model employed here, since the 
model explicitly contains all secular terms of the Hamilto- 
nian. Indeed, in accord with the reasoning presented above, 
throughout the duration of the integrations presented in the 
main text, the external perturber retains e^ = 0.5 while the 
eccentricities of the disk annuli never exceed e < 0.1 (see 
Figure 4). Furthermore, it is important to recall that the 
estimates of eccentricity excitation presented here should 
be viewed as an upper-limits, since dissipative forces of the 
gas will generally act to circularize the orbits on timescales 
much shorter than ~ 1 Myr. 

The inclination evolution of the considered disk is pre- 
sented in Figure 2. An important property to be noted from 
the inclination solution is that all plotted annuli never de- 
part significantly from each-other in phase-space. Instead, 
they always remain within i < few degrees of each-other. 
Physically, this means that the disk remains locally un- 
warped. An important consequence of this lack of warping 
is that standard theories of planetary formation"^^ and disk- 
driven migration^^ are directly applicable. In other words, 
because the disk's ascending node is precessing sufficiently 
slowly for the adiabatic condition^'^ to be satisfied, the pro- 
cesses of planetary formation and migration move forward 
as if the disk was completely isolated from the external 
perturber. 

While the eccentricities within the disk remain low 
throughout the duration of the integrations, the elliptic- 



Fig. 6. — Eccentricities of the disk annuli as functions of 
time corresponding to the Gaussian secular solutions pre- 
sented in Figure 5. See also Figure 4. 



ity of the stellar perturber's orbit does play an appreciable 
role in dictating the nodal recession frequency of the disk. 
Specifically, a higher e' corresponds to faster nodal reces- 
sion (in the subsequent sections, this will be shown to be 
a result of the non-linear secular terms contained in the 
Hamiltonian). To demonstrate this effect, as well as the 
sensitive dependence on inclination at near-orthogonal an- 
gles numerically, we performed an additional pair of nu- 
merical experiments (setting i' = 45 deg and i' = 75 
deg) with an identical setup to that described in the main 
text, except for a lower perturber eccentricity (e' = 0.1). 
The resulting inclination evolution is shown in Figure 5 
with the corresponding eccentricity evolution shown in Fig- 
ure 6. Note that although the mass and semi-major axis 
of the stellar perturber are held fixed, the changes in e' 
and i' make a significant difference in the nodal recession 
timescale compared to what is reported in Figure 2. Par- 
ticularly, the recession periods characteristic of the low ec- 
centricity (e^ = 0.1) setup are Tdisk — 1.4Myr (i' = 45 
deg) and Tdisk — 4.6Myr (i' = 75 deg), while that corre- 
sponding to the nominal (e' = 0.5) case are Tdisk — 0.9Myr 
{i' = 45 deg) and Tdisk - l-SMyr {i' = 70 deg). 

2. Laplace-Lagrange N-Ring Model 

As already mentioned above, the single-CPU implemen- 
tation of the Gaussian N-ring model is too computationally 
demanding to be a useful tool for mapping out parameter 
space. Fortunately, in the adiabatic regime where the disk 
remains coherent and unwarped, Laplace-Lagrange secular 
theory can be used to approximately reproduce the results 
obtained with the Gaussian model. As above, we model 
the proto-planetary disk as a conservative system of mas- 
sive concentric rings (initialized at i = 0), perturbed by an 
inclined distant massive ring which represents the stellar 
companion. 

The Laplace-Lagrange secular theory is qualitatively 
equivalent to the Gaussian secular model employed above, 
however the theory explicitly assumes that all eccentrici- 
ties and inclinations remain small, such that sin(i) ~ i. 
This assumption completely decouples the eccentricity and 
inclination dynamics. While neglecting the eccentricity- 
inclination coupling is justified for low-e^ perturbers, the 



resulting equations do not capture the dependence of the 
nodal recession rate on e' in case the latter is significant. 
Thus, for the purposes of this section, we shall restrict our 
discussion to low-e^ systems and focus on the analytical 
reproduction of the integrations shown in Figure 5. Sub- 
sequently, we shall introduce a leading-order correction for 
the perturber's eccentricity in the following section. 

In terms of Keplerian orbital elements, the scaled 
Laplace-Lagrange Hamiltonian of the j*^ ring reads: 
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JJ'j~^ Z^ BjkijikCos{Qj -Qk) (1) 

where as before, i is inclination, fl is the longitude of as- 
cending node and B's are interaction coefficients. While 
Keplerian orbital elements do not form a canonically con- 
jugated set, in terms of complex Poincare variables ^ = 
iexp(zr2), the Laplace-Lagrange equations of motion form 
an eigen-system" 
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where n denotes the mean motion, m is mass, a 
if (ttj < ak)]ak/aj if (a^ < aj), and ajk = djk if 
(ttj < ttk)] 1 if (ttk < aj). In contrast to planetary secular 
theory, in disk secular theory it is customary to soften the 
Laplace coefficients, ^3/25 by the disk aspect ratio h = H/a 
to account for its finite vertical thickness^^: 
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The closed-form solution to the equations of motion reads: 



0(0 = 5Z ^ok exp^(/fct + 5k) 



fc=i 



(5) 



where /'s are the eigen-frequencies and ^^'s are eigen-vectors 
of the B matrix, and (5's are phases which are set by initial 
conditions'^ . 

Because of the analytic nature of the solution, we are 
not as limited in N, allowing for a somewhat different rep- 
resentation of the disk. While, we take aout = 50AU as 
above, the inner edge is considerably closer to the host 
star: ain = 0.05AU. Furthermore, following ref. [41], the 
disk is now represented as a series oi N = 100 rings where 
the spacing is logarithmic with the semi-major axis ratio 
between neighboring annuli set to o; = 0.9325 through- 
out the disk. As before, the surface density profile of 
the disk is taken to be E oc r"-*^ and the disk mass is 

rridisk = lO-^Mstar. 

We make one trivial modification to the standard setup 
of the calculation: we reduce the coefficients in the B ma- 
trix corresponding to disk-external perturber interactions 
by a factor of cos(z^), treating only the disk's mid-plane pro- 
jected component of the force exerted by the external ring 
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Fig. 7. — Time evolution of the stellar-spin/disk orbit 
normal misalignment angle, -0. This solution was obtained 
using the Laplace-Lagrange secular theory and is aimed 
at reproducing Figure 5. The plotted of curves represent 
the annuli of the disk at a = lAU, a = lOAU and a = 
50AU. The nodal recession periods characteristic of this 
setup are Tdisk - 1.2Myr (i' = 45°) and Tdisk - 4Myr 
(i' = 75°). The agreement between the Laplace-Lagrange 
and Gaussian models is satisfactory, especially provided the 
level of approximation inherent to the Laplace-Lagrange 
theory (e.g. truncation of the disturbing function, e^ = 0, 
etc) and a somewhat different representation of the disk. 



as dynamically important. This introduces the l/cos(2^) 
dependence into the forced nodal recession period of the 
disk approximately observed in the Gaussian secular model. 

Figure 7 shows a reproduction of the examples shown 
above with Laplace-Lagrange secular theory. Aside from a 
mild quantitative difference between the observed nodal re- 
cession periods {Tj^^j^ ~ 1.2 Myr as compared to T^^^^^ — 
1.4 Myr for i' = 45 deg and T^^^^ — 4 Myr as compared 
to T^^j^^^ ^ 4.6 Myr for i' = 75 deg) the agreement be- 
tween the models is satisfactory, especially provided the 
level of approximation inherent to the Laplace-Lagrange 
theory and a somewhat different representation of the disk. 

At first glance, the applicability of classical Laplace- 
Lagrange secular theory in the context of the problem at 
hand may be surprising, given that the inclination between 
the disk and the external perturber is not small. However, 
this discrepancy is only apparent in the initial reference 
frame of the disk (i.e. the frame in which we measure ip). 
In a reference frame that is centered on the host star and 
is normal to the stellar orbit angular momentum vector (as 
shown in Figure 1), the disk's inclination, i^ , remains nearly 
constant in time. Rather, the quantity that changes is the 
disk's ascending node on which the Laplace-Lagrange the- 
ory places no restriction. It is further important to note 
that within the (self-gravitating) disk, the small-angle ap- 
proximation poses no problem since under adiabatic per- 
turbations considered here, the mutual inclinations among 
the neighboring annuli remain small at all times. Con- 
sequently, once the effect of the large mutual inclination 
between the disk and the stellar orbit is taken into account 
by projecting the force exerted by the external perturber 
onto the disk's mid-plane, Laplace-Lagrange secular the- 
ory provides an adequate approximation to the adiabatic 



dynamical evolution of the system. 

In summary, although a self-consistent account for 
non-secular perturbations as well as dissipative, non- 
gravitational effects within the disk will modify the solu- 
tion quantitatively, the qualitative behavior of the system 
is well captured within the context of the secular models 
used here. 

3. An Approximate Closed-Form Solution 

As already demonstrated with the Gaussian and 
Laplace-Lagrange models above, protoplanetary disks re- 
main largely unwarped thanks to their own self-gravity. 
We can take advantage of this fact to derive a simple, 
closed-form analytical solution for the secular evolution 
of continuous rigid disks under external perturbations. 
Recall that the evolution of the disk can simply be inter- 
preted as nodal recession at constant inclination. If all 
annuli of the disk are fixed to its mid-plane, the disk's 
state vector can be represented as a single point in polar 
coordinates with the radius r = \^\ = i^ and the polar 
angle 6 = arctan(Im[^]/Re[^]) = Q. Accordingly, a single 
recession period of the disk will be represented as a circle 
or radius r = i' centered on on the origin. 

The rate of nodal recession can be calculated using 
Laplace-Lagrange theory as above. In particular, the disk's 
recession rate is given by the orbital angular momentum 
weighted average of the forced recession rates of the disk 
annuli (i.e. the diagonal components of the B matrix). 
Upon letting N — )► oo, we can turn the sum over the rings 
that represent the disk into an integral of surface density: 
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Similarly, the recession rate of the perturbing star's node 
reads: 
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dt 
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The assumptions made up to this point generally imply 
that the recession rate of the perturbing body will be much 
slower than that of the disk itself and can be interpreted 
as the slow recession of the coordinate system in which the 
disk's state vector is measured. As a result, the period for 
cyclic excitation of the misalignment angle, ip, is given by 
Tdisk = |27r/(d (Qdzsk) /dt + dn'/dt)\. 

In the framework of the polar coordinates described 
above, the misalignment angle, -0, is represented by the 
phase-space distance between the disk's initial condition 
and its state vector. As the disk's node recesses, ip will 
trace the rim of the disk's circular phase-space orbit and 
its time evolution will therefore be described by a cycloid. 
Specifically, the evolution of ip is given by the following 
parametric equation: 

{t,^P} = {rd^sk{l - sin(7))/27r,i^(l - cos(7))} (8) 

where 7 is a parameter that governs the number of cycles 
through which the system has rotated (a single cycle cor- 
responds to 7 = 27r). We recalculated the evolution of tp 
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Fig. 8. — Time evolution of the stellar-spin/disk orbit 
normal misalignment angle, ip. This solution was obtained 
using the approximate cycloid model, based on the Laplace- 
Lagrange secular theory and is aimed at reproducing Figure 
5, and by extension. Figure 7 as well. The nodal recession 
periods characteristic of this setup are Tdisk — 1.3Myr (i' = 
45°) and Td^sk - 3.6Myr {i^ = 75°). 



for the example shown in Figures (5) and (7) with this sim- 
plified formalism, choosing the disk parameters as above. 
The calculated cycloids are presented in Figure 8. The 
observed nodal recession periods are quantitatively simi- 
lar for low inclination to those obtained with the Gaussian 
model (T^Jf. ^ 1.3 Myr as compared to T^^j^^^ ^ 1.4 
Myr for i' = 45 deg and T^^J^ ~ 3.6 Myr as compared to 
Y^Go^^ss ^ 4 g i^yj. fQj. ■/ ^ ^5 jgg^ j^ -g noteworthy that 
the Tdisk ex: cos(i) dependence is explicitly built into the 
equations. Furthermore, note that although the formalism 
is based on the Laplace-Lagrange theory, the characteristic 
periods of the cycloids differ somewhat from the Laplace- 
Lagrange N-ring solution due to the continuous represen- 
tation of the disk. 

In our discussion of the analytical models based on the 
Laplace-Lagrange theory so far, we have neglected the ef- 
fects of the coupling between the disk's inclination and the 
perturber's eccentricity. As already stated above, this as- 
sumption is only satisfied when e' is small (a setup not nec- 
essarily characteristic of star-forming environments). Con- 
sequently, here we will extend the model developed above 
to incorporate a correction for the perturber's eccentricity. 
To leading order, the coupling between the disk's inclina- 
tion and the perturber's eccentricity manifests itself as a 
fourth order secular effect. Accordingly, we write the new 
Hamiltonian of the disk as: 
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where the newly introduced interaction coefficient reads 
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Fig. 9. — Time evolution of the stellar-spin/disk orbit 
normal misalignment angle, ip. This solution was obtained 
using the approximate cycloid model, based on the Laplace- 
Lagrange secular theory, and corrected for the perturber's 
eccentricity, e^ The nodal recession periods characteristic 
of this setup are Tdisk — IMyr (i' = 45°) and Tdisk — 2Myr 
(i' = 70°), similar to those reported in Figure 2. 
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and k' refers to the index of the ring that represents the 
perturber. Generally, the fourth order secular Hamiltonian 
does not yield equations of motion that have analytical so- 
lutions. However, in the special case where most of the 
angular momentum of the system is contained in the or- 
bit of the stellar perturber (as is the case here), the per- 
turbers orbital elements remain nearly constant throughout 
the disk's lifetime. Consequently, e' can be treated as a pa- 
rameter rather than a variable, rendering the newly added 
coupling, an effectively second order term. 

Further simplifications can be made by taking advan- 
tage of the fact that ajy <C 1. Specifically, we expand 
the Laplace coefficient as a hypergeometric series, leading 



to the approximations h. 
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terms, we find that Bjj and Cjy take on similar forms. 
Consequently, after some rearrangement, we can rewrite 
the diagonal coefficients of the B matrix as 



^'33 



B, 



(11) 



and effectively reduce the Hamiltonian (9) to an expression 
that is formally identical to that of the Hamiltonian (1). 
This modification has a simple translation to the cycloid 
model introduced in this section: we can account for the 
perturber eccentricity by simply enhancing d (fldisk) /dt as 
given in equation (6) by a factor of (1 + 3e^^/2). 

Using the newly modified cycloid model, we can now 
aim to reproduce the results shown in Figure (2) of the 
main text quantitatively by explicitly setting e' = 0.5. The 
results are shown in Figure (9). The observed nodal reces- 
sion periods match those obtained with the Gaussian model 



disk 

Myr for i' = 45 deg and T^^J^ 



quite well: {T^^'^j^ ~ 1 Myr as compared to T^^j^^^ ^ 0.9 



n^Gaus 
disk 



1.8 Myr for i' 



2 Myr as compared to 
70 deg). Collectively, this anal- 



ysis suggests that significant eccentricity of the perturber 
can diminish the timescale for excitation of spin-orbit mis- 
alignment by a factor of ~ 2. 

4. Adiabatic Trailing of the Host Star 

The angular momentum of a host-star of mass Mstar, 
radius Rstar, and rotation rate uu can be represented as an 
orbiting ring with semi-major axis 



V 9PGMstar 



1/3 



and mass 



1/3 
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where ^2 is the apsidal motion constant, and / is the stel- 
lar moment of inertia. Upon making this approximation, 
we can interpret the forced nodal recession period of this 
ring (stellar bulge) as the characteristic angular momen- 
tum coupling timescale between the disk and the stellar 
spin axis, Tstar — 27t/B. As in the previous section, we 
consider a continuous disk by making A/" — )► co and take 
advantage of the fact that a <C ain, again leading to the 
approximation b^L — ^^- ^^^ resulting expression reads: 



'star '■ 



I G Mstar r«-* 
V a Ja^^ 



Ms 



rdr 
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We consider a typical pre-main-sequence star with M = 
IM0, Rstar = '2Rq, k2 = 0.014, / = 0.08, surrounded by 
the proto-planetary disk with an inner edge, ain, at the 
stellar co-rotation radius and outer edge at aout = 50 AU, 
comprising rridisk = 10~^ Mstar with a E oc r~^ surface 
density profile as above. This setup yields Tstar ~ 10 Myr 
and Tstar ~ 0.3 Myr for slow and fast rotators respectively. 
These estimates suggest that the spin axes of slow rotators 
will generally not adiabatically follow, bur rather slowly 
precess around their disk's angular momentum vectors. If 
Tdisk <^ Tstar, then the host-star's spin-axis will effectively 
remain stationary (as assumed in Figure 2 and Figure 5). 
However, if Tdisk ~ Ttar, significant excitation of mutual 
inclination between the disk and the host star can still take 
place, but the time-evolution of ip will be more complicated 
than the cyclic patterns presented in this work. 
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